2021 Fasting Cancer Project

Set options and working directory

Load libraries

library(Seurat)
library(cowplot)
library(umap)
library(dplyr)
library(Matrix)
library(readxl)
library(openxlsx)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(sctransform)
library(knitr)
library(kableExtra)
#library(biomaRt)
library(DESeq2)
library(escape)
library(dittoSeq)
library(GSEABase)
library(scater)
library(ComplexHeatmap)

Import marker sets

gene.lists <- read_xlsx("Munoz_Yilmaz_CellCycle_signatures.xlsx")
gene.lists.names <- colnames(gene.lists)
GOI.lists <- c()
for (i in gene.lists.names) {
  tmpList <- gene.lists %>% dplyr::select(all_of(i))
  tmpList <- tmpList[!is.na(tmpList)]
  GOI.lists[[i]] <- tmpList
}

Load the Cell Ranger Matrix Data and create the base Seurat object.*

the initial processing was done with r 3.6.1 with Seurat 3.2.0 - the UMAP comes out slightly differently in r 4.0.3 with Seurat 3.2.3*

#al.dat <- Read10X("200218Yil_data/al/filtered_feature_bc_matrix/")
#f.dat <- Read10X("200218Yil_data/f/filtered_feature_bc_matrix/")
#rf.dat <- Read10X("200218Yil_data/rf/filtered_feature_bc_matrix/")
#rf.rapa.dat <- Read10X("200218Yil_data/rf.rapa/filtered_feature_bc_matrix/")

#al <- CreateSeuratObject(counts = al.dat, project = "al", min.cells = 3, min.features = 200)
#f <- CreateSeuratObject(counts = f.dat, project = "f", min.cells = 3, min.features = 200)
#rf <- CreateSeuratObject(counts = rf.dat, project = "rf", min.cells = 3, min.features = 200)
#rf.rapa <- CreateSeuratObject(counts = rf.rapa.dat, project = "rf.rapa", min.cells = 3, min.features = 200)

#al[["percent.mt"]] <- PercentageFeatureSet(al, pattern = "^mt-")
#f[["percent.mt"]] <- PercentageFeatureSet(f, pattern = "^mt-")
#rf[["percent.mt"]] <- PercentageFeatureSet(rf, pattern = "^mt-")
#rf.rapa[["percent.mt"]] <- PercentageFeatureSet(rf.rapa, pattern = "^mt-")

#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(rf, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#VlnPlot(al, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

Merge samples to single Seurat Object

#merged_Seu <- merge(al, c(f,rf,rf.rapa), project = "Diet")

#merged_Seu <- NormalizeData(merged_Seu, normalization.method = "LogNormalize", scale.factor = 10000)
#merged_Seu <- FindVariableFeatures(merged_Seu, selection.method = "vst", nfeatures = 2000)
#merged_Seu <- ScaleData(merged_Seu)
#merged_Seu <- RunPCA(merged_Seu, features = VariableFeatures(object = merged_Seu))
#merged_Seu <- RunUMAP(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- RunTSNE(merged_Seu, reduction="pca",dims=1:30)
#merged_Seu <- FindNeighbors(merged_Seu, dims = 1:30, verbose = FALSE)
#merged_Seu <- FindClusters(merged_Seu, verbose = FALSE)

#DimPlot(merged_Seu,reduction="umap",group.by="orig.ident",label=TRUE,repel=FALSE)

#FeaturePlot(merged_Seu,reduction="umap",features="mt-Cytb",min.cutoff=0,max.cutoff=4,cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="percent.mt",cols=c("grey","red"))
#FeaturePlot(merged_Seu,reduction="umap",features="nFeature_RNA",cols=c("grey","red"))

Save/Load seurat object

#save(merged_Seu, file="merged.Robj")
#load("merged.Robj")

Dataset Integration

#di <- SplitObject(merged_Seu, split.by = "orig.ident")

#for (i in 1:length(di)) {
#  di[[i]] <- NormalizeData(di[[i]], verbose = FALSE)
#  di[[i]] <- FindVariableFeatures(di[[i]], selection.method = "vst", nfeatures = 2000,
#                                      verbose = FALSE)
#}

#dat.anchors <- FindIntegrationAnchors(object.list=di,dims=1:30)
#integrated <- IntegrateData(anchorset=dat.anchors,dim=1:30)

#DefaultAssay(integrated)<-"integrated"
#integrated <- ScaleData(integrated)
#integrated <- RunPCA(integrated,npcs=30)
#integrated <- RunUMAP(integrated,reduction="pca",dims=1:30)
#integrated <- RunTSNE(integrated,reduction="pca",dims=1:30)
#integrated <- FindNeighbors(integrated, dims = 1:30, verbose = FALSE)
#integrated <- FindClusters(integrated, verbose = FALSE)

Load integrated data and create UMAP from original integrated run*

NOTE: loading final object to avoid recalculating ssGSEA data

#load("integrated_orig.Robj")
load("../repo_data/integrated_091120.Robj")
#load("../repo_data/integrated_final.Robj")

Test plots

DefaultAssay(integrated)<-"integrated"
DimPlot(integrated,reduction="umap",split.by="orig.ident",group.by="orig.ident")

DimPlot(integrated,reduction="umap",group.by="orig.ident")

DimPlot(integrated,reduction="umap",group.by="integrated_snn_res.0.8", label=TRUE)

## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.

Figure 3B - Add cluster/treatment metadata columns and plot labeled UMAP

integrated@meta.data$treat_clust <- paste0(integrated@meta.data$orig.ident,integrated@meta.data$integrated_snn_res.0.8)
integrated@meta.data$clust_treat <- paste0(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
integrated@meta.data$celltype <- integrated@meta.data$integrated_snn_res.0.8
Idents(object = integrated) <- integrated$celltype

integrated <- RenameIdents(object = integrated,  '16' = '16_Tuft')
integrated <- RenameIdents(object = integrated,  '11' = '11_EC')
integrated <- RenameIdents(object = integrated,  '13' = '13_EEC')
integrated <- RenameIdents(object = integrated,  '2' = '2_Stem')
integrated <- RenameIdents(object = integrated,  '5' = '5_Stem')
integrated <- RenameIdents(object = integrated,  '10' = '10_Stem')
integrated <- RenameIdents(object = integrated,  '14' = '14_Paneth')
integrated <- RenameIdents(object = integrated,  '8' = '8_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '9' = '9_Goblet')
integrated <- RenameIdents(object = integrated,  '1' = '1_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '4' = '4_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '15' = '15_Secretory_Progenitor')
integrated <- RenameIdents(object = integrated,  '3' = '3_EC_Progenitor')
integrated <- RenameIdents(object = integrated,  '6' = '6_EC_Progenitor')
integrated <- RenameIdents(object = integrated,  '0' = '0_Early_TA')
integrated <- RenameIdents(object = integrated,  '7' = '7_Early_TA')
integrated <- RenameIdents(object = integrated,  '12' = '12_Unknown')


integrated[["clust_celltype"]] <- Idents(object = integrated)

Fig3b UMAP Plot with color blind safe colors

Idents(object = integrated) <- integrated$clust_celltype

celltype_colors <- c("2_Stem"="#117733",
                     "5_Stem"="#999933",
                     "10_Stem"="#009E73",
                     "0_Early_TA"="#E69F00",
                     "7_Early_TA"="#D55E00",
                     "6_EC_Progenitor"="#0072B2",
                     "3_EC_Progenitor"="#56B4E9",
                     "11_EC"="#88CCEE",
                     "13_EEC"="#6699CC",
                     "1_Secretory_Progenitor"="#661100",
                     "4_Secretory_Progenitor"="#882255",
                     "8_Secretory_Progenitor"="#CC6677",
                     "15_Secretory_Progenitor"="#AA4499",
                     "14_Paneth"="#332288",
                     "16_Tuft"="#000000",
                     "9_Goblet"="#F0E442",
                     "12_Unknown"="#999999")  

dp.cb <- DimPlot(integrated,reduction="umap", cols=celltype_colors, label=TRUE, repel=TRUE, pt.size=2, label.size=6) + NoLegend()
dp.cb

#pdf('Fig3b.cbSafe.pdf',width=14, height=10)
#dp.cb
#dev.off()

Figure 3C - LGR5 Vln plot with color blind safe colors

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
plotOrder <- c("5","2","10","0","1","3","4","6","7","8","9","11","12","13","14","15","16")

vln_colors <- c("2"="#117733",
                "5"="#999933",
                "10"="#009E73",
                "0"="#E69F00",
                "1"="#661100",
                "3"="#56B4E9",
                "4"="#882255",
                "6"="#0072B2",
                "7"="#D55E00",
                "8"="#CC6677",
                "9"="#F0E442",
                "11"="#88CCEE",
                "12"="#999999",
                "13"="#6699CC",
                "14"="#332288",
                "15"="#AA4499",
                "16"="#000000")  

Idents(integrated) <- factor(Idents(integrated), levels= plotOrder)
vl.cb <- VlnPlot(integrated, cols=vln_colors, idents = , features = c("Lgr5"), pt.size = 0.5, slot="data")
vl.cb

#pdf('Fig3c_cbSafe.pdf',width=14, height=8)
#vl.cb
#dev.off()

Figure S3,B - Table of cell counts in each integrated data cluster and sample

p.integrated <- table(integrated@meta.data$integrated_snn_res.0.8,integrated@meta.data$orig.ident)
round(prop.table(p.integrated,2),3)

##     
##         al     f    rf rf.rapa
##   0  0.169 0.191 0.135   0.178
##   1  0.096 0.101 0.135   0.119
##   2  0.082 0.116 0.104   0.102
##   3  0.099 0.091 0.098   0.101
##   4  0.125 0.078 0.082   0.061
##   5  0.076 0.080 0.084   0.072
##   6  0.076 0.074 0.065   0.077
##   7  0.056 0.040 0.096   0.055
##   8  0.078 0.053 0.034   0.037
##   9  0.028 0.046 0.043   0.060
##   10 0.034 0.032 0.044   0.055
##   11 0.026 0.056 0.033   0.029
##   12 0.039 0.028 0.024   0.032
##   13 0.005 0.006 0.009   0.009
##   14 0.005 0.005 0.006   0.007
##   15 0.002 0.002 0.005   0.005
##   16 0.004 0.002 0.004   0.001

Figure 5D - heatmap

DefaultAssay(integrated)<- "integrated"
Idents(object = integrated) <- integrated$integrated_snn_res.0.8
dd <- subset(integrated, idents = c("2", "5", "10"), downsample=100)

topvst <- head(VariableFeatures(dd), 500)
mat <- as.matrix(dd@assays$integrated@scale.data) #as.matrix(subset_dd@assays$integrated@scale.data)
mat <- mat[topvst,]

genes <- c(GOI.lists$Biton_S1_ISC.I, GOI.lists$Biton_S1_ISC.II, GOI.lists$Biton_S1_ISC.III)
labels <- c(rep('Biton ISC I', length(GOI.lists$Biton_S1_ISC.I)), 
            rep('Biton ISC II', length(GOI.lists$Biton_S1_ISC.II)), 
            rep('Biton ISC III', length(GOI.lists$Biton_S1_ISC.III)))

idxs <- which(genes %in% rownames(mat))
genes <- genes[idxs]
labels <- labels[idxs]
mat <- mat[genes,]

ht <- Heatmap(mat, column_names_side = 'top', row_names_gp = gpar(fontsize = 9), column_names_gp = gpar(fontsize = 12),
              column_title = '', name = 'scaled data', column_labels = rep('', ncol(mat)),
              column_split = factor(as.character(dd$integrated_snn_res.0.8), levels = c('5', '2', '10')), 
              row_split = labels, row_order = genes, #column_order = sort(colnames(mat)),
              cluster_column_slices = FALSE,
              top_annotation = HeatmapAnnotation(cluster = as.character(dd$integrated_snn_res.0.8)))

draw(ht)

#pdf('Fig3D.pdf',width=12, height=10)
#draw(ht)
#dev.off()

GSEA dot plots

data <- as.data.frame(read_excel('cluster_2_5_10stem_gsea.xlsx', sheet = 'Sheet2'))
data

##                NAME                Comparison        NES    FDR q-val
## 1    BITON_S1_ISC.I         f5_v_al5.noNA.NES  1.3215407 0.2682069200
## 2   BITON_S1_ISC.II         f5_v_al5.noNA.NES  2.0391748 0.0000000000
## 3  BITON_S1_ISC.III         f5_v_al5.noNA.NES -1.1991149 0.3889467000
## 4    BITON_S1_ISC.I        rf5_v_al5.noNA.NES  1.7629898 0.0019750001
## 5   BITON_S1_ISC.II        rf5_v_al5.noNA.NES  2.2566988 0.0000000000
## 6  BITON_S1_ISC.III        rf5_v_al5.noNA.NES  0.8156711 0.9565237000
## 7    BITON_S1_ISC.I   rf.rapa5_v_al5.noNA.NES  0.9629948 0.9024616000
## 8   BITON_S1_ISC.II   rf.rapa5_v_al5.noNA.NES  2.0106223 0.0000000000
## 9  BITON_S1_ISC.III   rf.rapa5_v_al5.noNA.NES  0.7934456 0.9918505500
## 10   BITON_S1_ISC.I         f2_v_al2.noNA.NES  1.3257033 0.2230794400
## 11  BITON_S1_ISC.II         f2_v_al2.noNA.NES  2.0813794 0.0000000000
## 12 BITON_S1_ISC.III         f2_v_al2.noNA.NES -0.8689138 0.6766603600
## 13   BITON_S1_ISC.I        rf2_v_al2.noNA.NES  1.2952416 0.3336595600
## 14  BITON_S1_ISC.II        rf2_v_al2.noNA.NES  2.0858900 0.0000000000
## 15 BITON_S1_ISC.III        rf2_v_al2.noNA.NES  0.9813662 0.9102961000
## 16   BITON_S1_ISC.I   rf.rapa2_v_al2.noNA.NES -1.9235135 0.0047468350
## 17  BITON_S1_ISC.II   rf.rapa2_v_al2.noNA.NES  1.8073893 0.0003529412
## 18 BITON_S1_ISC.III   rf.rapa2_v_al2.noNA.NES  1.2005840 0.4414399000
## 19   BITON_S1_ISC.I       f10_v_al10.noNA.NES  1.0625067 0.6857741500
## 20  BITON_S1_ISC.II       f10_v_al10.noNA.NES  2.1364486 0.0000000000
## 21 BITON_S1_ISC.III       f10_v_al10.noNA.NES  0.7806379 1.0000000000
## 22   BITON_S1_ISC.I      rf10_v_al10.noNA.NES  0.8879838 0.8637682000
## 23  BITON_S1_ISC.II      rf10_v_al10.noNA.NES  2.0704706 0.0000000000
## 24 BITON_S1_ISC.III      rf10_v_al10.noNA.NES  0.6868493 0.9426303500
## 25   BITON_S1_ISC.I rf.rapa10_v_al10.noNA.NES  0.9626449 0.9918727000
## 26  BITON_S1_ISC.II rf.rapa10_v_al10.noNA.NES  2.0112374 0.0000000000
## 27 BITON_S1_ISC.III rf.rapa10_v_al10.noNA.NES  0.9090225 0.9317172000

data$FDR <- data$`FDR q-val` + 0.001
data$Comparison <- gsub('\\.noNA\\.NES','',data$Comparison)

comparisons <- c('f5_v_al5', 'rf5_v_al5', 'rf.rapa5_v_al5')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl5 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 5 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl5

#pdf('Fig5F.pdf',width=4, height=4)
#cl5
#dev.off()

comparisons <- c('f2_v_al2', 'rf2_v_al2', 'rf.rapa2_v_al2')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl2 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 2 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl2

#pdf('FigS3Da.pdf',width=4, height=4)
#cl2
#dev.off()

comparisons <- c('f10_v_al10', 'rf10_v_al10', 'rf.rapa10_v_al10')
sub_data <- data[which(data$Comparison %in% comparisons),]
cl10 <- ggplot(data=sub_data, aes(y=NAME, x=factor(Comparison, levels = comparisons), size=-log10(FDR), color=NES)) + 
  geom_point() + 
  scale_color_gradient2(midpoint=0, low="blue", mid="white",
                        high="red", space ="Lab", limits=c(-3,3))+
  scale_size_continuous(range=c(1,6))+
  ggtitle('Cluster 10 GSEA') + theme_classic() + 
  theme(legend.position="right", axis.text.x = element_text(angle = 90)) + ylab('Gene Set') + xlab('Comparison')
cl10

#pdf('FigS3Db.pdf',width=4, height=4)
#cl10
#dev.off()

#there is a glitch in this plot, cl2 loses x axis, legend order is different

figS3D <- ggarrange(cl2,cl10, ncol=2, nrow=1)
figS3D

#pdf('FigS3D_incorrect.pdf',width=8, height=4)
#figS3D
#dev.off()

Extended data Figure S4A - Feature Plots - color blind safe

DefaultAssay(integrated)<-"RNA"
colorScheme <- c("#C1BEBF","#fe6100")

fp.Muc2 <- FeaturePlot(integrated,reduction="umap",features="Muc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Muc2

fp.Tff3 <- FeaturePlot(integrated,reduction="umap",features="Tff3",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Tff3

fp.Lyz1 <- FeaturePlot(integrated,reduction="umap",features="Lyz1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lyz1

fp.Defa30 <- FeaturePlot(integrated,reduction="umap",features="Defa30",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Defa30

fp.Chga <- FeaturePlot(integrated,reduction="umap",features="Chga",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Chga

fp.Reg3a <- FeaturePlot(integrated,reduction="umap",features="Reg3a",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Reg3a

fp.Alpi <- FeaturePlot(integrated,reduction="umap",features="Alpi",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Alpi

fp.Atoh1 <- FeaturePlot(integrated,reduction="umap",features="Atoh1",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Atoh1

fp.Lgr5 <- FeaturePlot(integrated,reduction="umap",features="Lgr5",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Lgr5

fp.Smoc2 <- FeaturePlot(integrated,reduction="umap",features="Smoc2",cols=colorScheme, slot = "data",label=TRUE, repel=TRUE,label.size=4)
fp.Smoc2

figS4A <- ggarrange(fp.Muc2,fp.Tff3,fp.Lyz1,fp.Defa30,fp.Chga,fp.Reg3a,fp.Alpi,fp.Atoh1,fp.Lgr5,fp.Smoc2, legend = "right", ncol = 5, nrow = 2)

figS4A

#pdf('FigS4A_cbSafe.pdf',width=20, height=10)
#figS4A
#dev.off()

Figure 5E - cell cycle data - color blind safe

DefaultAssay(integrated)<-"RNA"
gene.names <- rownames(integrated@assays$RNA@data)
Idents(object = integrated) <- integrated$seurat_clusters
stem.subset <- subset(integrated, idents = c("2","5","10"))
levels(stem.subset) <- c("5","2","10")

vln_colors <- c("2"="#117733",
                "5"="#999933",
                "10"="#009E73")
  
s.InData <- intersect(gene.names,GOI.lists$mm.s)
g2m.InData <- intersect(gene.names,GOI.lists$mm.g2m)
stem.subset[["percent.mm.s"]] <- PercentageFeatureSet(stem.subset, features = s.InData)
stem.subset[["percent.mm.g2m"]] <- PercentageFeatureSet(stem.subset, features = g2m.InData)
mm.s <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.s",pt.size = 0.3,slot = "data")
mm.s

mm.g2m <- VlnPlot(stem.subset, cols = vln_colors, features="percent.mm.g2m",pt.size = 0.3,slot = "data")
mm.g2m

figS4b <- ggarrange(mm.s,mm.g2m, legend = FALSE, ncol=2, nrow=1)
figS4b

pdf('figS4b_cbSafe.pdf',width=8, height=4)
figS4b
dev.off()

## png 
##   2

integrated[["percent.mm.s"]] <- PercentageFeatureSet(integrated, features = s.InData)
integrated[["percent.mm.g2m"]] <- PercentageFeatureSet(integrated, features = g2m.InData)

mm.g2m.Flag <- if_else(integrated@meta.data$percent.mm.g2m >= 0.3, "Yes", "No")
integrated@meta.data$mm.g2m.Flag <- mm.g2m.Flag

p <- table(integrated$clust_treat,integrated$mm.g2m.Flag)
p.g2m.summary <- round(prop.table(p,2),3)

mm.s.Flag <- if_else(integrated@meta.data$percent.mm.s >= 0.2, "Yes", "No")
integrated@meta.data$mm.s.Flag <- mm.s.Flag

p <- table(integrated$clust_treat,integrated$mm.s.Flag)
p.s.summary <- round(prop.table(p,1),3)

escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters

This code is no longer working due to R and package updates but resulting data is stored in the seurat object escape 1.0.0 is probably required but is no longer available. the version in this R image is 1.0.1 escape ssGSEA - run ssGSEA to quantify expression of the BitonI gene set clusters

#egs <- GeneSet(GOI.lists$Biton_S1_ISC.I, setName="BitonI")
#ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
#integrated@meta.data$BitonI_ssGSEA <- ES$BitonI

escape ssGSEA - run ssGSEA to quantify expression of the BitonII gene set clusters

# egs <- GeneSet(GOI.lists$Biton_S1_ISC.II, setName="BitonII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonII_ssGSEA <- ES$BitonII

escape ssGSEA - run ssGSEA to quantify expression of the BitonIII gene set clusters

# egs <- GeneSet(GOI.lists$Biton_S1_ISC.III, setName="BitonIII")
# ES <- enrichIt(obj = integrated, gene.sets = egs, groups = 1000, cores = 4)
# integrated@meta.data$BitonIII_ssGSEA <- ES$BitonIII

Extended data 5 plots - colorblind safe

Figure 5G - Biton 1 in cluster5 Figure 5G_II - Biton II in cluster2 Figure 5G_III - Biton III in cluster10 Figure 5H - Pdgfa in cluster5 Figure S3E - Gkn3 in cluster5*

Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust

vln_colors <- c("al5"="#e69f00",
                "f5"="#56b4f9",
                "rf5"="#117733",
                "rf.rapa5"="#d55e00")

plot.order <- c("al5","f5","rf5","rf.rapa5")
Idents(object = clust5.subset) <- factor(Idents(object = clust5.subset), levels = plot.order)
vln.BitonI.ssGSEA <- VlnPlot(clust5.subset, cols = vln_colors, features="BitonI_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-I_ssGSEA")

vln.BitonI.ssGSEA

#pdf('Fig5G.pdf',width=4, height=4)
#vln.BitonI.ssGSEA
#dev.off()

vln.Pdgfa <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Pdgfa"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Pdgfa

pdf('FigS4g_cbSafe.pdf',width=4, height=4)
vln.Pdgfa
dev.off()

## png 
##   2

vln.Gkn3 <- VlnPlot(clust5.subset, cols = vln_colors, features = c("Gkn3"), slot= "data",pt.size = 0.4) + NoLegend()
vln.Gkn3

pdf('FigS4h_cbSafe.pdf',width=4, height=4)
vln.Gkn3
dev.off()

## png 
##   2

clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust

vln_colors <- c("al2"="#e69f00",
                "f2"="#56b4f9",
                "rf2"="#117733",
                "rf.rapa2"="#d55e00")

plot.order <- c("al2","f2","rf2","rf.rapa2")
Idents(object = clust2.subset) <- factor(Idents(object = clust2.subset), levels = plot.order)
vln.BitonII.ssGSEA <- VlnPlot(clust2.subset, cols = vln_colors, features="BitonII_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-II_ssGSEA")
vln.BitonII.ssGSEA

#pdf('Fig5G_II.pdf',width=4, height=4)
#vln.BitonII.ssGSEA
#dev.off()

clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust

vln_colors <- c("al10"="#e69f00",
                "f10"="#56b4f9",
                "rf10"="#117733",
                "rf.rapa10"="#d55e00")

plot.order <- c("al10","f10","rf10","rf.rapa10")
Idents(object = clust10.subset) <- factor(Idents(object = clust10.subset), levels = plot.order)
vln.BitonIII.ssGSEA <- VlnPlot(clust10.subset, cols = vln_colors, features="BitonIII_ssGSEA",slot = "data",pt.size = 0.4) + 
  stat_summary(fun = median, geom='point', size = 15, colour = "black", shape = 95) + NoLegend() + ylab("Enrichment Score") + ggtitle("Biton ISC-III_ssGSEA")
vln.BitonIII.ssGSEA

#pdf('Fig5G_III.pdf',width=4, height=4)
#vln.BitonIII.ssGSEA
#dev.off()

Biton plots assembled

bitonPlot <- ggarrange(vln.BitonI.ssGSEA,vln.BitonII.ssGSEA,vln.BitonIII.ssGSEA,nrow=1,ncol=3)
bitonPlot 

#pdf('FigS4f_cbSafe.pdf',width=10, height=5)
#bitonPlot
#dev.off()

Fig 5H style plots for Oat

Figure 5H style plots for Oat, Oct, Ass Asl in cluster5

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.subset <- subset(integrated, idents = c("5"))
Idents(object = clust5.subset) <- clust5.subset$treat_clust
levels(clust5.subset) <- c("al5","f5","rf5","rf.rapa5")
clust5.subset$treat_clust <- factor(x = clust5.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5"))
VlnPlot(clust5.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.

vln.Oat.5 <- VlnPlot(clust5.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5

#pdf('Fig5H_Oat.5.pdf',width=4, height=4)
#vln.Oat.5
#dev.off()

clust2.subset <- subset(integrated, idents = c("2"))
Idents(object = clust2.subset) <- clust2.subset$treat_clust
levels(clust2.subset) <- c("al2","f2","rf2","rf.rapa2")
clust2.subset$treat_clust <- factor(x = clust2.subset$treat_clust, levels = c("al2","f2","rf2","rf.rapa2"))
VlnPlot(clust2.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

vln.Oat.2 <- VlnPlot(clust2.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.2

#pdf('Fig5H_Oat.2.pdf',width=4, height=4)
#vln.Oat.2
#dev.off()

clust10.subset <- subset(integrated, idents = c("10"))
Idents(object = clust10.subset) <- clust10.subset$treat_clust
levels(clust10.subset) <- c("al10","f10","rf10","rf.rapa10")
clust10.subset$treat_clust <- factor(x = clust10.subset$treat_clust, levels = c("al10","f10","rf10","rf.rapa10"))
VlnPlot(clust10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

#vln.Oat.10 <- VlnPlot(clust10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
#vln.Oat.10

#pdf('Fig5H_Oat.10.pdf',width=4, height=4)
#vln.Oat.10
#dev.off()

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$seurat_clusters
clust5.2.10.subset <- subset(integrated, idents = c("5","2","10"))
Idents(object = clust5.2.10.subset) <- clust5.2.10.subset$treat_clust
levels(clust5.2.10.subset) <- c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10")
clust5.2.10.subset$treat_clust <- factor(x = clust5.2.10.subset$treat_clust, levels = c("al5","f5","rf5","rf.rapa5","al2","f2","rf2","rf.rapa2","al10","f10","rf10","rf.rapa10"))
VlnPlot(clust5.2.10.subset, features = c("Oat"), split.by = "treat_clust", group.by = "treat_clust", slot= "data",pt.size = 0.5)

vln.Oat.5.2.10 <- VlnPlot(clust5.2.10.subset, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.5.2.10

#pdf('Fig5H_Oat.5.2.10.pdf',width=12, height=4)
#vln.Oat.5.2.10
#dev.off()

Figure 5H style plots but without cluster for Oat*

DefaultAssay(integrated)<-"RNA"
Idents(object = integrated) <- integrated$orig.ident
levels(integrated) <- c("al","f","rf","rf.rapa")
integrated$orig.ident <- factor(x = integrated$orig.ident, levels = c("al","f","rf","rf.rapa"))
VlnPlot(integrated, features = c("Oat"), split.by = "orig.ident", group.by = "orig.ident", slot= "data",pt.size = 0.5)

vln.Oat.noclust <- VlnPlot(integrated, features = c("Oat"), slot= "data",pt.size = 0.1)
vln.Oat.noclust

#pdf('Fig5H_Oat.noclust.pdf',width=4, height=4)
#vln.Oat.noclust
#dev.off()

Figure S3E - Differential Expression Testing to prepare .rnk files*

# DefaultAssay(integrated)<-"RNA"
# Idents(object = integrated) <- integrated$treat_clust
# levels(x = integrated)
# ##
# #adjust the id1 and id2 variables to set up different tests
# #could use a loop for this
# ##
# id1 <- "al5"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
# 
# id1 <- "al10"
# id2 <- "al2"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)
# 
# id1 <- "al10"
# id2 <- "al5"
# outData <- paste(id1,id2,"Markers",sep=".")
# outData <- FindMarkers(integrated, ident.1 = id1, ident.2 = id2,logfc.threshold = 0, assay="RNA")
# outData <- tibble::rownames_to_column(outData,"#MGISym")
# rnk.tmp <- outData %>% dplyr::select(c("#MGISym","avg_logFC"))
# #write.table(rnk.tmp, sep='\t',file=paste0(id1,"_v_",id2,".rnk"),col.names=TRUE, quote=FALSE, row.names=FALSE)

Save integrated object with ssGSEA data

#save(integrated, file="integrated_final.Robj")

sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] ComplexHeatmap_2.6.2        scater_1.18.6              
##  [3] SingleCellExperiment_1.12.0 GSEABase_1.52.1            
##  [5] graph_1.68.0                annotate_1.68.0            
##  [7] XML_3.99-0.5                AnnotationDbi_1.52.0       
##  [9] dittoSeq_1.2.6              escape_1.0.1               
## [11] DESeq2_1.30.1               SummarizedExperiment_1.20.0
## [13] Biobase_2.50.0              MatrixGenerics_1.2.1       
## [15] matrixStats_0.58.0          GenomicRanges_1.42.0       
## [17] GenomeInfoDb_1.26.7         IRanges_2.24.1             
## [19] S4Vectors_0.28.1            BiocGenerics_0.36.1        
## [21] kableExtra_1.3.1            knitr_1.31                 
## [23] sctransform_0.3.2           ggpubr_0.4.0               
## [25] ggrepel_0.9.1               ggplot2_3.3.3              
## [27] openxlsx_4.2.3              readxl_1.3.1               
## [29] Matrix_1.2-18               dplyr_1.0.4                
## [31] umap_0.2.7.0                cowplot_1.1.1              
## [33] Seurat_3.2.3               
## 
## loaded via a namespace (and not attached):
##   [1] reticulate_1.18           tidyselect_1.1.0         
##   [3] RSQLite_2.2.3             htmlwidgets_1.5.3        
##   [5] BiocParallel_1.24.1       Rtsne_0.15               
##   [7] munsell_0.5.0             codetools_0.2-16         
##   [9] ica_1.0-2                 future_1.21.0            
##  [11] miniUI_0.1.1.1            withr_3.0.0              
##  [13] colorspace_2.0-0          highr_0.8                
##  [15] rstudioapi_0.13           ROCR_1.0-11              
##  [17] ggsignif_0.6.0            tensor_1.5               
##  [19] listenv_0.8.0             labeling_0.4.2           
##  [21] GenomeInfoDbData_1.2.4    polyclip_1.10-0          
##  [23] farver_2.0.3              pheatmap_1.0.12          
##  [25] bit64_4.0.5               parallelly_1.23.0        
##  [27] vctrs_0.6.5               generics_0.1.0           
##  [29] xfun_0.21                 R6_2.5.1                 
##  [31] clue_0.3-58               ggbeeswarm_0.6.0         
##  [33] rsvd_1.0.3                msigdbr_7.2.1            
##  [35] locfit_1.5-9.4            bitops_1.0-6             
##  [37] spatstat.utils_2.0-0      cachem_1.0.3             
##  [39] DelayedArray_0.16.3       assertthat_0.2.1         
##  [41] promises_1.1.1            scales_1.1.1             
##  [43] beeswarm_0.2.3            gtable_0.3.0             
##  [45] beachmat_2.6.4            Cairo_1.5-12.2           
##  [47] globals_0.14.0            goftest_1.2-2            
##  [49] rlang_1.1.4               genefilter_1.72.1        
##  [51] GlobalOptions_0.1.2       splines_4.0.3            
##  [53] rstatix_0.6.0             lazyeval_0.2.2           
##  [55] broom_0.7.4               yaml_2.2.1               
##  [57] reshape2_1.4.4            abind_1.4-5              
##  [59] backports_1.2.1           httpuv_1.5.5             
##  [61] tools_4.0.3               ellipsis_0.3.2           
##  [63] RColorBrewer_1.1-2        ggridges_0.5.3           
##  [65] Rcpp_1.0.6                plyr_1.8.6               
##  [67] sparseMatrixStats_1.2.1   zlibbioc_1.36.0          
##  [69] purrr_1.0.2               RCurl_1.98-1.2           
##  [71] rpart_4.1-15              openssl_2.2.0            
##  [73] deldir_0.2-9              GetoptLong_1.0.5         
##  [75] viridis_0.5.1             pbapply_1.4-3            
##  [77] zoo_1.8-8                 haven_2.3.1              
##  [79] cluster_2.1.0             magrittr_2.0.3           
##  [81] data.table_1.13.6         RSpectra_0.16-0          
##  [83] scattermore_0.7           circlize_0.4.12          
##  [85] lmtest_0.9-38             RANN_2.6.1               
##  [87] fitdistrplus_1.1-3        GSVA_1.38.2              
##  [89] hms_1.0.0                 patchwork_1.1.1          
##  [91] mime_0.9                  evaluate_0.24.0          
##  [93] xtable_1.8-4              rio_0.5.16               
##  [95] shape_1.4.5               gridExtra_2.3            
##  [97] compiler_4.0.3            tibble_3.0.6             
##  [99] KernSmooth_2.23-17        crayon_1.4.1             
## [101] htmltools_0.5.8.1         mgcv_1.8-33              
## [103] later_1.1.0.1             tidyr_1.1.2              
## [105] geneplotter_1.68.0        DBI_1.1.1                
## [107] MASS_7.3-53               car_3.0-10               
## [109] cli_3.6.3                 igraph_1.2.6             
## [111] forcats_0.5.1             pkgconfig_2.0.3          
## [113] foreign_0.8-80            scuttle_1.0.4            
## [115] plotly_4.9.3              xml2_1.3.2               
## [117] vipor_0.4.5               webshot_0.5.2            
## [119] XVector_0.30.0            rvest_0.3.6              
## [121] stringr_1.4.0             digest_0.6.36            
## [123] RcppAnnoy_0.0.18          spatstat.data_2.0-0      
## [125] rmarkdown_2.6             cellranger_1.1.0         
## [127] leiden_0.3.7              edgeR_3.32.1             
## [129] uwot_0.1.10               DelayedMatrixStats_1.12.3
## [131] curl_5.2.1                shiny_1.6.0              
## [133] rjson_0.2.20              lifecycle_1.0.4          
## [135] nlme_3.1-149              jsonlite_1.8.8           
## [137] BiocNeighbors_1.8.2       carData_3.0-4            
## [139] limma_3.46.0              viridisLite_0.3.0        
## [141] askpass_1.1               pillar_1.4.7             
## [143] lattice_0.20-41           fastmap_1.2.0            
## [145] httr_1.4.2                survival_3.2-7           
## [147] glue_1.4.2                zip_2.1.1                
## [149] spatstat_1.64-1           png_0.1-7                
## [151] bit_4.0.4                 stringi_1.5.3            
## [153] blob_1.2.1                BiocSingular_1.6.0       
## [155] memoise_2.0.1             irlba_2.3.3              
## [157] future.apply_1.7.0

writeLines(capture.output(sessionInfo()), "2024_sessionInfo.txt")